Procedure for enhancing resolution of spectral information

ABSTRACT

A procedure for enhancing the resolution of spectral data, in which the Fourier self-deconvolution method FSD is used to produce from the input spectrum a data set in an interval 0-L t , the maximum entropy method MEM is used to compute prediction error filter coefficients {a k  } from this data, and by the linear prediction method LP, using coefficients {a k  } and data points 0 to M. Data are predicted in the interferogram I(x) beyond L t , whereby output spectrum maximum line narrowing with minimum distortion is achieved.

The present invention concerns a procedure performed entirely in the time domain, for increasing the resolution of spectral information, or for reducing the width of spectral lines, in FTIR, IR, NMR or other spectroscopy.

For instance when the quantities of gaseous impurities are monitored in industrial areas, factory halls, switching yards or other critical locations, or when air quality in general is studied, one has to be able to observe concentrations on the order of a few ppb. It is typical, in apparatus constructed for this purpose, to use chemical or semiconductor detectors which are sensitive to a single gas or a few gases only. Since infrared spectroscopy, at its best, enables the same sensitivity to be reached, it becomes possible to monitor simultaneously almost all measurable gases, and this being furthermore a comparatively rapid method of analysis, it is a noteworthy alternative in applications of this type.

However, the spectra obtained as results of measurement may be difficult of their interpretation. The gas which is being measured may contain many different constituents and various interference factors which impede the interpretation of the spectrum. The object of the present invention is to disclose a design which enables more accurate and more detailed interpretation of measured spectral information than before, and elimination of disturbances of the measurement.

Fourier Self-Deconvolution (FSD) and the Maximum Entropy Method (MEM) are both well-established procedures serving spectral line stabilizing (spectral information resolution enhancing). FSD is at present mainly employed in infrared .spectroscopy of solid or liquid state, while MEM has been applied in numerous other branches of science as well.

FSD constitutes a linear method of deconvolution which is exclusively based on multiplication of the measured information in the time domain. In the present application text the FT-IR convention is followed, referring to the interferogram domain as the time domain and the spectral domain, as the frequency domain. MEM, on the other hand, represents a non-linear fitting method maximizing the information entropy of the output data. It is typical of this method that minor changes at the input cause powerful distortions in the output, and the results may even depend on the algorithm used. Furthermore, MEM gives only one of a plurality of possible solutions, whereas FSD yields an exact solution containing noise. While MEM is a well-approved procedure today, many practitioners of this method add a warning that it generates ghost lines or fails to find real, small lines.

In the following we shall first compare the FSD and MEM procedures as regards their efficiency to narrow spectral lines.

We shall use the following pair of Fourier transformations in this text: ##EQU1## where E(ν) is the spectrum (i.e., intensity as a function of wave number), I(x) is the interferogram (i.e., intensity as a function of optical path difference), and { } and ⁻¹ { } are the Fourier transform and the inverse Fourier transform, respectively.

Fourier self-deconvolution (FSD) takes as input a spectrum E(ν), consisting of lines with assumed line shapes W(ν). These input line shapes are converted to a desired output line shape W_(o) (ν) by the time domain operation specified in Formula 2: ##EQU2## where W'(ν) is the deconvolved spectrum having the output line shape W_(o) (ν). If the desired line shape is narrower than the initial line shape, line narrowing has been achieved. The line shape, line narrowing has been achieved. The line narrowing factor K, loosely referred to as "resolution enhancement factor", is given by the formula: ##EQU3## where FWHH(W) and FWHH(W_(o) ) are the full width at half height of W(ν) and W_(o) (ν), respectively.

The main advantage of FSD is that it is a linear procedure, that is, the Fourier self-deconvolutions of E₁ (ν) and E₂ (ν):

    E.sub.1 (ν)+E.sub.2 (ν)=E.sub.1 (ν)+E.sub.2 (ν)(4)

It is well-known (notably in electronics) that an output signal S' (t) can be defined as the convolution of the input signal S(t) and the impulse response h(t) (of a circuit). Therefore, in the time domain, ##EQU4## Using the interferometric nomenclature, let us postulate that the input signal I(x)= ⁻¹ {E(ν)} in the time domain has passed through the filter defined by the impulse response h(t). Then, the output signal will be:

    Γ(x)=h(x)*I(x)                                       (6)

This can be expressed in discrete form as ##EQU5## where M specifies the length of the impulse response function and h_(k) in the convolution sum contains the sampling interval x. The error {e_(j) } between the input data {I_(j) } and output data {I'_(j) } is given by ##EQU6## The digital filter defined by the impulse response {a_(k) } k=8,1,2,. . .M is called the prediction-error filter because the output of this filter is the error vector {e_(j) }. Notice that the filter defined by {h_(k) } is different from {a_(k) }. Using Equations 6 and 7, the input spectrum will be given by

    E(ν)= {I(x)}= {Γ(x)+e(x)}= {h(x)*I(x)+e(x)}= {h(x)}.E(ν)+ (ν)(10)

where (ν) is now the error spectrum. From Eq. 10 the input spectrum E(ν) can then be derived as ##EQU7## In this formula, h(x) or a(x) can be computed (see above), but (ν) must be approximated in some way. From Eq. 11, with the denominator implicitly expressed by the use of Eq. 9, the power spectrum of E(ν) can be estimated as ##EQU8## In other words, σ_(M) ² is simply the expectation value of

the error power, or σ_(M) is the standard deviation of the error vector {e_(k) }. According to Equations 8 and 9, the input data {I_(j) } obey the following equation: ##EQU9## where N is the number of sampling points. Upon multiplying the formula 14 I^(*) _(j-n) (complex conjugate of I_(jn)), and taking the expectation value (i.e., the average) for both sides of this equation, one obtains for 0≦j-n≦N and 0≦j-k≦N: ##EQU10## Since e_(j) and I^(*) _(j-n) do not correlate if n=/0, we have ##EQU11## If the autocorrelation coefficient R_(n) is ##EQU12## Eq. 15 can be expressed in the following form: ##EQU13## Since a_(o) =1, Eq. 18 becomes ##EQU14## which is known as the Yule-Walker equation and can be expressed in matrix form as ##EQU15##

The maximum entropy spectral estimation in its basic form is accomplished as follows. First, the {a_(k) } coefficients are derived by use of Eq. 19 or 20 so as to minimize σ_(M) ². Then, when the power spectrum of E( ) is computed from Eq. 12, the information entropy S=-Constant, *Σ|E_(j) |ln|E_(j) | is maximized. Thus, minimizing σ_(M) ² in Eq. 20, will maximize the entropy S, provided the approximation in Eq. 12 is valid.

However, the solution to Eq. 20 poses some problems. A major difficulty is that the autocorrelation coefficients R_(k) are not well defined with the aid of Eq. 17, owing to the finite length of the data set. This problem has been eliminated by Burg's modification of MEM, in which the coefficients {a_(k) } k=l,2,. . .,M, as well as σ_(M) ⁻¹² and R_(k), are recursively computed from the previous values {a_(k) } k=l,2,. . . ,M-1, σ_(M-1) ² and R_(k), so that the error power, representing the average of the forward and backward predictions, is minimized.

The linear prediction (LP) method is based on Eq. 14. Since the computed coefficients {a_(k) } minimize σ_(M) and the error vector {e_(j) }, we can approximate e_(j=I) _(j) -I'_(j) =0 or I_(j=I') _(j). Now Eq. 14 can be rewritten as ##EQU16## If {I_(j) } j=0,1,2,. . .,M and the coefficients {a_(k) } j=l,2,3,. . .,M are known, the data set may be extended beyond I_(M) by recursively computing new data I_(M+1) as follows: ##EQU17##

Since I_(j) =I'_(j), we later use a denotation where the original data set is the input signal, and the calculated data set is the output signal. Thus, in Linear Prediction, we assume that the first M+1 output data points are equal to the M+1 input data points. Beyond the M^(th) data point there exist only the (predicted) output data points.

In the following, the technique of prior art and the procedure of the invention shall be described in detail, referring to the attached drawings, wherein FIG. 1, top row, shows a test spectrum of three Lorentzian lines without noise, Column A is the FSD spectrum with a Bessel output line shape when K=2, 3, 4 and 5, and Column B is the spectrum computed by Burg's MEM method with prediction error filters M=6, 8, 12 and 15 and the data points numbering 50,

FIG. 2 is the same presentation as in FIG. 1, but with a test spectrum having signal/noise ratio 1000,

FIG. 3 illustrates the noise in FSD when K=5, and the distortion in MEM when M=15, the top row being the bottom row of FIG. 2,

FIG. 4 presents a comparison of the FSD and MEM procedures in a case of extremely noisy spectra (S/N=3),

FIG. 5 is a comparison of FSD and FSD+MEM in a case with favourable signal/noise ratio (S/N=1000): Column A is the same as in FIG. 2, while Column B is the spectrum obtained by MEM from the FSD spectrum, with sinc line shape,

FIG. 6 a diagrammatic presentation of the LOMEP method along the route indicated by solid arrows, the dotted line arrows correlating the data in the time domain with the corresponding frequency domain presentation,

FIG. 7 is a comparison of the FSD and LOMEP spectra in the case of poor signal/noise ratio (S/N=10), and

FIG. 8 displays the advantages of the trifunction method LOMEP over FSD, compared in the case of a spectrum with poor S/N ratio, the top spectrum being composed of Lorentzian lines at 200, 600 and 610 cm⁻¹ of equal intensities and widths.

In the following, the Fourier self-deconvolution method shall be compared with the maximum entropy method, and for test spectrum three Lorentzian lines are used which have equal widths and intensities and where one line is isolated and the other two are overlapping.

In FIG. 1 is presented the comparison of the FSD and MEM methods in a case in which there is no noise in the test spectrum. Column A in FIG. 1 depicts the FSD spectrum, which has Bessel output line shape, W_(o) (ν) in FIG. 2, and where the line narrowing coefficient K=2, 3, 4 and 5. Column B depicts the amplitude spectrum computed by Burg's MEM method with prediction error filters M=6, 8, 12 and 15. The amplitude spectrum was derived from formula 12 and taken as the square root of the power spectrum. For a given number of input data points the number M is the sole adjustable parameter in MEM. The top row of each column in FIG. 1 shows the input spectrum. Column A illustrates the FSD method of line narrowing, which is applied to the test data. The line shapes of the deconvoluted spectrum are precisely Bessel. On the other hand, the MEM method, illustrated in Column B, yields an output spectrum with various distortions depending on the number M of the prediction error filter. It is fair, in this case, to point out that the FSD method uses the information of the line shape, whereas MEM does not.

In FIG. 2 the above comparison of FSD and MEM is carried on to a case in which random noise has been added to the input spectrum so that the signal/noise ratio is 1000. Once again, the FSD method results in a spectrum having the same degree of line narrowing as in FIG. 1A, but with increased noise. It is important to note that the noise level (r ms) is constant over the spectrum and that there is no correlation between lines and noise, i.e., the noise on a line is the same as the noise elsewhere. This means that the peak heights never deviate from their theoretical values by more than the noise amplitude. The situation is different in the spectrum of Column B in FIG. 2. The output spectrum obtained by MEM has no noise, only distortion of the lines. In other words, in the regions without lines the MEM spectrum is noise-free, yet strong distortions are present close to the true lines, which are impossible to estimate. Also, the degree of line narrowing is not readily adjustable in MEM by changing M.

FIG. 3 illustrates the noise/distortion in FSD and in MEM. The lower row in this figure displays in FSD, when K=5, the accrued noise and in MEM, when M=15, the incurred total distortion. The upper row is the bottom row of FIG. 2.

A comparison of the FSD and MEM methods in the case of an extremely noisy spectrum (S/N not more than 3) is shown in FIG. 4. When FIG. 3B is compared with FIG. 2B, it is obviously evident that the nature of the MEM spectrum's distortion has changed, but there is not yet any evidence of white noise. On the other hand, the spectrum obtained by the FSD method (FIG. 4A) reveals a considerable deterioration in the signal-to-noise ratio.

The question arises as to why MEM does not yield better results than what is shown in FIG. 1B. The reason is that the original MEM method of Burg rests on the basic assumption that the time series must be stationary, i.e., not decaying. In other words, the MEM approach postulates that the signal S(t) consists of a sum (see the formula) in which the prediction error filter acts like an electronic filter. This, however, is not the case in solid or liquid state IR or NMR spectra. The time domain signals have a roughly exponential decay, owing to their Lorentzian line shape. In addition, Eq. 12 contains an approximation, namely |E(ν)|² =σ_(M) ², which implies that even when the {a_(k) } coefficients are correct, |E(ν)|² can still be distorted.

The maximum entropy method is best applied to spectra consisting only of sinusoidal signals with constant amplitudes in the time domain. Accordingly, the best line shape for MEM would be a sine function (i.e., the Fourier transform of a boxcar function). This fact inspired the Applicant to use FSD to change the actual line shape to a sinc function, i.e., W_(o) (ν) now becomes a sinc function in Eq. 2.

We shall start with a combination of FSD and MEM. FIG. 5 shows spectra obtained with Burg's modified MEM, but derived from an FSD spectrum with sinc line shapes. In Column A in FIG. 5 are depicted the output spectra obtained by FSD (see FIG. 2A), in order to facilitate comparison. In the FSD+MEM method, the number M is the number of non-zero data in the corresponding interferogram. This obviously results in a greatly improved MEM output, which is not unexpected since by using FSD data for input to MEM the line shape information is now taken into account, which should improve the output data. However, there are still distortions to the intensity and shape of the lines, and the line narrowing process is not well controlled with MEM.

The object of the present invention is to eliminate the above-described drawbacks of methods of prior art. Regarding the features which are characteristic of the invention, reference is made to the claims section.

The invention discloses a novel procedure which combines FSD with some features of MEM and with the corollary thereof, linear prediction (LP). This combination of FSD, MEM and LP is called the LOMEP method (Line shape Optimized Maximum Entropy linear Prediction). The essence of the procedure of the invention is utilizing of the strong points of each of these methods in order to achieve maximum line narrowing with minimum distortion.

The distortions in the MEM spectrum are due to Eq. 12. Therefore, instead of using the formula 12, line narrowing is accomplished by applying linear prediction to extend the signal. In the time domain there exists a point x=L_(t) (t=truncated) where the signal becomes roughly equal to the noise. We apply MEM to the "good" data between 0 and L_(t) only to derive the correlation coefficients {a_(k) } needed in linear prediction, formula 22.

Thus, in the trifunction method LOMEP the FSD function is first used to produce a set of data in the interval from 0 to L_(t). The second function MEM is then used to generate the coefficients {a_(k) } from these data. Finally, the third function LP uses the coefficients {a_(k) } towards predicting the data in the interferogram beyond L_(t). FIG. 6 schematically illustrates the basic concept of the LOMEP method. The first step uses FSD to generate an interferogram I(x), consistent with a sinc function as output line shape. This interferogram is then used as input to MEM, from which the prediction error filter coefficients {a_(k) } are computed. In the present example the number M of the prediction error filter has been selected to equal the number of input data points, i e , L_(t) =M x; because the signals are now cosine waves with constant amplitudes, there is a strong correlation between all data points.

The next step is prediction of the extra data points in the interferogram I(x) beyond the initial point L_(t), by using the coefficients {a_(k) } and the data points from 0 to M. This is accomplished by means of linear prediction, given in the formula 22. Theoretically, the Fourier transform of this expanded interferogram has sinc line shape. Therefore, prior to computing the final spectrum, the expanded interferogram should be apodized in order to give a more useful line shape in the final spectrum. The ultimate K factor obtained with LOMEP depends on L_(f) =N_(f) x (where N_(f) +1 is now the total number of data points) and on the desired output line shape.

FIG. 7B illustrates the results that have been obtained with the LOMEP method, applied to a spectrum with S/N=10. As shown in Column A, the maximum line narrowing realistically achievable in this case with FSD is K=2. On the other hand it is possible with the trifunction method LOMEP to attain K values of 5 and higher.

Furthermore, the S/N ratio of the LOMEP output spectrum is not dependent on the K value: it depends only on L_(t), which in its turn is a function of the initial S/N. The M value used in FIG. 7B was 60, and the S/N ratio is roughly the same for all K values, the reason being that the noise only comes from the data in the interval from 0 to L_(t) in the interferogram. There is no noise in the predicted portion x>L_(t). It should also be noted that the position, shape and intensity of the lines are only affected by the noise, as in the case of the FSD spectrum. Consequently, this method retains all the advantages of the linear methods. Again, since only the data between 0 and L_(t) are used, LOMEP does not alter the intensity, nor the position, of the input lines. The data predicted beyond L_(t), on the other hand, only changes (i.e., reduces) the line widths if the coefficients {a_(k) } are correct. Compare for instance the spectrum of FIG. 7B with the spectrum presented in FIGS. 5A, 2B and 1B, all of which have different S/N ratios.

In the last example, FIG. 8 illustrates a case in which two lines in the neighbourhood of 600 cm⁻¹ lie so close to each other that FSD is not successful in their complete separation, whereas LOMEP leads to a well-resolved doublet with correct intensities (peak heights).

The critical factor in optimation of LOMEP is the value of L_(t), which is a function of the S/N ratio in the original spectrum. A remaining question is, how many extra data points should be computed by formula 22. In theory this should not be a critical factor; in practice, however, noise and/or an erroneous FSD line shape will set a certain limit to this number. 

I claim:
 1. A method of examining a sample to determine the presence or absence of a component therein comprising:producing a spectral image of a sample in the form of spectral data, spectroscopically producing spectral data from a sample representing the composition of said sample, producing from said spectral data an interferogram of said sample between 0 and L_(f) ; producing by the Fourier self-deconvolution method FSD a data set from the obtained interferogram in an interval from 0 to L_(t) where L_(t) <L_(f), said data set comprising data points O to M where M is the number of data points in said data set; computing from said data set by the maximum entropy method MEM prediction error filter coefficients {a_(k) }, where k=0,1,2,. . .M; predicting data beyond L_(t) in the interferogram using said prediction error filter coefficients {a_(k) }, and the data points 0 to M in said data set to achieve maximum line narrowing with minimum distortion in the resulting spectrum, and examining lines in said spectrum to detect presence or absence of a component in the composition of said sample, whereby the narrowing of the lines and the minimization of distortion facilitates the examination and detection of said component.
 2. A method according to claim 1, comprising applying said FSD to form said interferogram I(x) whose output line shape is consistent with a sinc function.
 3. A method according to claim 1, comprising, prior to computing the final spectrum apodizing the expanded interferogram I(x) to alter the line shape in the final spectrum.
 4. A method according to claim 1, comprising determining the extra data points in the linear prediction by the formula: ##EQU18##
 5. A method according to claim 1, comprising establishing the point x=L_(t) in the time domain to corresponds to stays in the point where the signal is substantially equal to the noise.
 6. A method according to claim 1, comprising rendering said interferogram substantially free of noise.
 7. A method according to claim 1, comprising providing the spectral data with a signal/noise ratio of 1000, 100, 10 or about
 3. 8. A method according to claim 1, wherein said sample is solid, liquid or gaseous. 